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We study the dynamics of Bose-Einstein condensates flowing in optical lattices on the basis of 
quantum field theory. For such a system, a Bose-Einstein condensate shows a unstable behavior 
which is called the dynamical instability. The unstable system is characterized by the appearance 
of modes with complex eigenvalues. Expanding the field operator in terms of excitation modes 
including complex ones, we attempt to diagonalize the unperturbative Hamiltonian and to find its 
eigenstates. It turns out that although the unperturbed Hamiltonian is not diagonalizable in the 
conventional bosonic representation the appropriate choice of physical states leads to a consistent 
formulation. Then we analyze the dynamics of the system in the regime of the linear response 
theory. Its numerical results are consitent with as those given by the discrete nonlinear Schrodinger 
equation. 

PACS numbers: 03.75.Lm,03.70.Kk,11.10.-z 



I. INTRODUCTION 

The Bose-Einstein condensates (BECs) of trapped 
atoms, first realized in 1995 [H, 0, S|, are ideal systems 
for studying quantum many-body phenomena. This is 
because the systems are dilute, weakly interacting ones 
and we can easily control the configuration of the trap 
and even the strength of atomic interaction. Recently 
much attention has been focused on the BEC in an opti- 
cal lattice, which provides us with rich phenomena such 
as superfluid-insulator quantum phase transition Q and 
Bloch oscillations and also enables us to observe 
directly the quantum fluctuation as quantum depletion 
from the condensate Q . 

Fluctuations and excitation spectra of a BEC are de- 
termined from the Bogoliubov-de Gennes (BdG) equa- 
tion, which is obtained by the linearization of the time 
dependent Gross-Pitaevskii (TDGP) equation. Numeri- 
cal calculation of the excitation spectra of a BEC in an 
optical lattice was performed in Ref. The analyti- 
cal studies under the tight-binding approximation [3, 
gave the theoretical result in a good agreement with the 
experiment 

In this paper, we focus on the BECs flowing in optical 
lattices. For such a system, two types of instability are 
known: Landau and dynamical instabilities. The Lan- 
dau instability is an energetic instability which is caused 
by negative energies of the quasi-particle The dy- 
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namical instability is that of condensates caused by the 
complex eigenvalues of the BdG equations [ll|, [13] ■ The 
Landau instability (energetic instability) is not observed 
at low temperature [l^, because it requires the thermal 
cloud which brings a dissipative mechanism and drives 
the condensate toward a lower energy state. On the other 
hand, the dynamical instability can occur even at zero 
temperature. 

The emergence of complex eigenvalues in the BdG 
equations is not restricted to the superfluidity of BECs in 
optical lattices, but is found widely, such as in the BECs 
with highly quantized vortices 
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solitons or in the multi-component BECs 
In all the cases, a naive scenario of the dynamical in- 
stability is as follow: When the complex eigenvalues ap- 
pear, the small deviations from the static condensation 
have the time-dependence with the complex frequency 
and will grow exponentially, and consequently the initial 
configuration of the static condensate is destroyed. The 
dynamical instability of the superfluidity flowing in an 
optical lattice was observed at the experiments and the 
decay rates of the condensate are well reproduced by the 
theory based on the TDGP equation [13, '22'|. In the case 
of a highly quantized vortex of BECs, the imaginary part 
of the complex eigenvalue is considered to be related to 
the experimentally observed lifetime of the vortex [Ts'l . It 
is also mentioned that the numerical result of the TDGP 
is in a good agreement with the experimental data [23| . 

The c-number theories seem to describe the dynami- 
cal instability well. But there are quantum fluctuations 
even at zero-temperature. So a consistent treatment of 
the BEC in full quantum field theory (QFT) is needed 
for describing the unstable behavior with quantum fluc- 
tuations. Because of the presence of trapping potentials, 
this QFT becomes one for finite volume systems in which 
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the spatially translational invariance is lost and all the 
energy levels are discrete. The Bose-Einstein conden- 
sation is considered as a spontaneous breakdown of a 
global phase symmetry, so the Nambu-Goldstone (NG) 




mode inevitably appears 24 1. This zero-mode requires 

2E 113]: closely related to its 



us subtle treatments [2 

observable effects [2^ [29[ and several theoretical matters 
such as the inequivalent vacua, Ward-Takahashi relations 
and Hugenholtz-Pines theorem [s^, [U [H, . 

It is desirable to include complex modes into the QFT 
in addition to the zero-mode. We have attempted to for- 
mulate QFT in the presence of complex-energy modes 
in the case of highly quantized vortices in the previous 
paper [3^ . In this paper we present a consistent formu- 
lation of QFT with complex modes for the BECs flowing 
in an optical lattice . These formulations in full quan- 
tum theory are expected to go beyond the conventional 
c-number approach of the mean-field approximation and 
the BdG equations. 

Needless to say, QFT in nonequilibrium situations is 
needed in wide area, for example, in particle physics, cos- 
mology and condensed matter physics, and its construc- 
tion is still a challenging subject. The QFT description 
of unstable behaviors, presented here, may become im- 
portant for some nonequilibrium systems. 

This paper is organized as follows. In Sec. II, the model 
action and Hamiltonian are given. We assume that the 
solution of the Gross-Pitaevskii (GP) equation satisfies 
the Bloch condition. In Sec. Ill, we expand the field 
operator in terms of the adequate complete orthonor- 
mal set and make the tight-binding approximation in the 
unperturbed Hamiltonian. In Sec. IV, solving the equa- 
tion of time evolution for the operators, we find eigen- 
modes, some of whose eigenvalues may be complex. We 
show that the operators in the complex mode sector are 
not subject to the usual bosonic commutation relations 
and that the unperturbed Hamiltonian can not be diag- 
onalizcd in the usual bosonic representation. One can 
check that the canonical commutation relations for the 
field operators which are the fundamental requirements 
in QFT are kept under the tight-binding approximation 
and in the presence of complex modes. Next, we ob- 
tain the eigenstates of complex modes and their proper- 
ties. In Appendix A, we summarize the eigenfunctions of 
the Bogoliubov-de Gennes (BdG) equation, in connection 
with the discussion in Sec. IV. In Sec. V, we introduce 
the conditions for physical states which provide us with a 
consistent QFT description of the unstable behavior. In 
Sec. VI, applying Kubo's linear response theory (LRT) to 
the physical states introduced in Sec V, we calculate the 
density response of the system against the external per- 
turbation, whose lengthy expression is given in Appendix 
B. Our numerical results of the density response repro- 
duce those of the discrete nonlinear Schrodinger equa- 
tion (DNSE), which is obtained by applying the tight- 
binding approximation to the TDGP equation [ssl. 13^. 
Section VII is devoted to summary. 



II. MODEL HAMILTONIAN OF QUANTUM 
FIELD THEORY FOR EEC IN OPTICAL 
LATTICE 

We consider the trapped BEG of neutral atoms in an 
optical lattice. In QFT, the model action which describes 
the system is given by 

S = j dtd^x^^\x){T - K -Vopt+ ii)-^{x) 

-^^\x)^\x)^{x)^{x)^ , (1) 



with x — (x,t). Here we use the following notations: 



K . -^V^ 
2m 



t—x,y,z 



Vopt{x) = Yj ^oCOS^ ( 



(2) 
(3) 
(4) 



The parameters /i, m and g represent the chemical po- 
tential, the mass of the neutral atom and the coupling 
constant of atomic interaction, respectively. The opti- 
cal potential whose strength is denoted by Vq has the 
lattice spacing di {i = x, y, z) in each direction, that is, 
Vopt{x + G) = Vopt{x), where G is a Bravais lattice vec- 
tor, G = {lxdx,lydy,lzdz) with integers U (i = x,y,z). 
We consider the situation without the harmonic potential 
for atoms. 

The action ([1]) is invariant under the global phase 
transformation ^i{x) e'^'^{x) and '^'^{x) —> e~'^\I't(a;), 
where 9 is an arbitrary real constant. When a BEG is cre- 
ated, the global phase symmetry is spontaneously broken. 
In the terminology of the operator formalism (canonical 
formalism) for quantum field theory (QFT), the Heisen- 
berg field '^{x) is then divided into a c-number part v{x) 
and an operator one 4'{x), 



^{x) = v{x) + (i){x) . 



(5) 



The c-nunibcr field v{x), called the order parameter, is 
assumed to be time- independent throughout this paper. 
The order parameter is defined as an expectation value 
of the Heisenberg field with respect to the vacuum 



or equivalently 



(r2|*(a;)|l]) = v(x), 



{n\^{x)\n) = 0. 



(6) 



(7) 



Let us introduce an additional symmetry breaking 
term [13 

AS^ j dtd^xel[v*{x)-^[x)+v{x)¥{x)\ , (8) 
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and the total action is defined by 
S, = S + AS. 



(9) 



Here e is an infinitesimal dimensionless parameter and 
e represents a typical energy scale of the system. The 
total action is not invariant under the global phase 
transformation of "^{x), due to the additional breaking 
term. 

The above method corresponds to the Bogoliubov's 
quasi-average, known in the treatment of spontaneously 
broken systems [13] ■ The breaking term specifies the "di- 
rection" of the symmetry breakdown, so a corresponding 
vacuum is selected among many degenerate ones. The 
singularity associated with the Nambu-Goldstone (NG) 
mode (23j is regularized [131 ■ At the final stage of calcu- 
lation, the limit of e — > is taken, and original symmetry 
is restored. 

We move to the interaction representation of the 
canonical formalism. The canonical commutation rela- 
tions (CCRs) for the field operators are given as follows: 

[$ix,t),^^{x',t)]^S{x-x'), (10) 
[${x,t),${x',t)] = [$^{x,t),^^{x',t)] =0. (11) 

The total Hamiltonian of the system is divided into the 
two terms. 



int : 



(12) 



where the unperturbative Hamiltonian Ho and the inter- 
action one ifint are given as 



(fx 



(13) 



d-'x 



V* (K + V - fi + g\v\^ ~ ee) 
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ee) V 



(14) 



respectively. The renormalization counter terms are sup- 
pressed. 

The condition ([7]) at the tree level (zero-loop level) 
derives the following classical equation for v{x): 



{K + Vopt+g\v{x)\^~H~ee)v{x)=0. 



(15) 



This equation is nothing but the Gross-Pitaevskii (GP) 
equation 5^ at the limit of e — > 0. The quantity |w(a;)p 
is identified with the density of condensed particles n{x) 
as n{x) — \v{x)\'^. The total condensate particle number 
A^c is given by A^c = / x\v{x)\'^. 



A. Bloch condition for Gross-Pitaevskii equation 

In this paper we assume that the density of condensate 
particles n{x) = |i;(a;)p has the same periodicity as the 



lattice as n{x + G) = n{x), although it may not be pe- 
riodic strictly. Under the assumption of this periodicity, 
the solution of the GP equation ([TS]) is expressed by a 
function satisfying the Bloch condition [38| : 



Vk{x) 



Ah-x 



Vk{x) 



Vkix + G) = Vkix) . 



(17) 



Hereafter we specify the solution of the GP equation by 
the vector k as Vk{x). Here U {i = x,y,z) ia an integer 
and Li represents a length of the system in each direction. 

In this paper, we are interested in the solution of the 
GP equation with the condensate flow in an optical lat- 
tice. Only when Vkix) is complex, the condensate can 
have a flow. The velocity of the condensate flow is given 
by 



Vk = —k -I \Sk{x) , 

m m 

where Sk{x) is the phase oivk{x). 



(18) 



III. TIGHT-BINDING APPROXIMATION IN 
OPTICAL LATTICE 

Here we discuss BECs in a cubic lattice for simplicity 
{d ~ dx ~ dy = dz). Consider the following eigenequa- 
tion: 



[K + Kpt + gW? - M - = €^frA^) 



(19) 



The density of the condensate particles nk{x) ~ \vk{x)\'^ 
has the periodicity of the lattice spacing d, so that one 
may require f'j^^ [x) to satisfy the Bloch condition as 



e^'^-/(J(a;) 



(20) 
(21) 



where q and n are a Bloch wave vector and a band index, 
respectively. Equation (|19p has the zero-energy solution 
which is proportional to Vk{x) as 



Vk{x) 



(22) 



with e^^2 = 0- The completeness and orthonormal condi- 
tions are 



d^xf^^/{x)f'j:^{x) = 5„,„V.' (23) 



hq 



5{x 



(24) 



n.q 



We consider the situation in which the strength of the 
lattice potential Vg is so large that the wave packet of 
the condensate wavefunction is localized in each site of 
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the lattice. The wave packet locahzed in each site is rep- 
resented by a Wannier function, whose set {wl^^ (x)} is 

related to the set {fj^ix)} by the unitary transforma- 
tion, 

4:'(^)=^s-^E^"^''"*/ff(^)' (25) 

where Ig is the total number of lattice sites and Xi is 
the position of the i-th site, where i = {ix,iy,iz) with 
integers ij {j = x,y,z). The Bloch wave vector q runs 
over the Brillouin zone. The completeness and orthonor- 
mal condition of the Wannier functions is guaranteed by 
Eqs. (ESI) and (121, 

jd^xw^^P*{x)w'^^^{x) = (26) 
E"'ii^*(^')4-H^) = 6{x-x'). (27) 

If the density of condensate particle nk{x) has the reflec- 
tion symmetry {nk{—x) = nk{x)), we can easily derive 
the following relation 

ft2;ix)^fSix). (28) 

Then it is shown from Eq. (^5]) that the Wannier func- 
tions become real. But, in our discussions below, we do 
not assume the reflection symmetry. 

The field operator is expanded in terms of the set 

<^(a;,t)=^4':)(t)4':)(a;), (29) 

71, i 

I 




Ho 



with the notations of 

j(nin2) / j3 (ni)* /r^ , ta \ ("2) 

^(«in2n3"4) ^ /d^ ("2) («3)* (n4) 

(36) 



where the operator a^" satisfies the usual bosonic 
commutation relations, '^Lj- = SijSnn' and 

[4"^4?] = [4':^^4f ^] = 0. The solution of the GP 
equation Vk{x) can be rewritten as 



Vkix) = nlY^e^'-^^wi'^x), (30) 

i 

Nc , . 

nc = ^. (31) 



Here we use Eq. (|22p and the inverse transformation of 

Eq. m. 



/ff(^) = ^s-^E^''-^^4:^(^)' (32) 



following from the formula. 



-rE^'^''"'^'"=Vq- (33) 



Substituting Eq. into Eq. and Eqs. and (1^ 
into Eq. (|13p , and making use of Eq. (|26p , one can rewrite 
the GP equation and the unperturbed Hamiltonian as 



(34) 



(35) 



Suppose that the system is so cold that we may restrict 
excitations only in the first band. Then the completeness 
condition holds in the first band. 



Y,wUx')wki{x)^b(x-x'), (37) 



.(n2) 
^fci2 



[(_ T(«in2) _ c. . c \ '("i)t^( 



ife-(xi,+Xi,) ~(ni)~(n2) 
■ "fcii "fei2 
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or equivalently from Eq. (|32p . 



(38) 



Here and hereafter the band index (n = 1) is omitted. 

By taking the most relevant terms for the kinetic terms 
Jiji2 and the on-site terms for interaction Ui-^i^i^i^ (tight- 
binding Umit), Eqs. (|34p and ((35)) are reduced to 

z-^i = J2 {Jie'^''^ + J!e-'''^'^)-Un, + ee, (39) 

l—x.y,z 



^0 - - ^ Jija\i 



o-kj 



^^{z~ fj. + 2U nc)al-aki 



UHc 



(40) 

where < i,j > represents the sum over the nearest 
neighbours and the following notations are introduced: 
Ji = Jii-\-in Ji — Ji-\-iii: z — J-a and U = Una- 
Here 1; is the unit vector along the Z-direction, e.g. 
Ix — (1,0,0). We consider the Fourier representation 
of the operators Uki and a^^^ 

aki = /7*^Cfcqe^'^-=% (41) 

4 = is^T.^Ue'"'''^ (42) 



where the operators Ckq and c^^ satisfy the bosonic com- 
mutation relations 



[cfcq,Cfeq'] = [cig,cl ,] =0. 



(43) 
(44) 



Substituting the Fourier representation of the opera- 
tors aki, CLki Sind the chemical potential (j39p into the 
unperturbed Hamiltonian we reduce it to 



-^0 - ^^kqci 



Uric 



kq'^kq 



EH 



kq^kq + Ck-qCkq + cl._^cl,- 



(45) 

Xkq- = Yl MJi\sn^(Qi + kid+^)sin(^)+e-e, 



l—x^y,z 



(46) 



where we write J; — | J/|e'®' and replace q with k + q. 
Hereafter our notation is simplified as A^q = Xkk+q ~^ 
Afcq. Each component of qi runs only over — 7r/d < 



(Ji < n/d {i = x,y,z) in the summation, because Ckq{t) 
and Afcq have the periodicity as Ckq+K{i) = Ckq{t) and 
Afcq+ic — Xkq where K is the reciprocal vector. 

The expansion of the quantum field operator Eq. 
can be rewritten with help of Eqs. p2|) and (|4ip as 



(47) 



Yckq{t)fkq{x) 



under the tight-binding approximation. 

IV. HAMILTONIAN WITH COMPLEX 
EIGENVALUES AND ITS EIGENSTATES 

Let us consider the equation of time evolution for 
Cfcq(t) and cl_-{t) as 

ifl^Ckqit) ^ Ckqit),Ho =TkqCkqit), (48) 

where we have introduced the doublet notation 

-'"'^(IS))- '''' 

The 2 X 2-matrix Tkq is given by 



Tkq — 



_ ( i>^kq + Unc 



-UUc 



Uric 
i^k-q + Uric) 



(50) 



To solve Eq. (|48p. we attempt to diagonalize Tkq- Let us 
set up the eigenequation 

TkqXkq = fT^kqXkq , (51) 

which gives two independent eigenvalues fno^ and dou- 
blet eigenvectors x^^ . The eigenvalues ftw^^'' are found 
to be 



where 



A 



kq 



^(•^fcq - ^k-q) , 
A|q + 2C/ncAfcq, 

^kq + ^k-q 



(52) 
(53) 
(54) 

(55) 



Note that these eigenvalues become complex if the con- 
dition, 



A^q + 2UncKkq < , 



(56) 

is satisfied. In Fig. [U we plot the stability phase diagram 
for different values ofncU/\J\. 

For a one-dimensional system, this condition is rewrit- 
ten simply as 



cos(fcd + 6) s cos(fcc? -I- 8) sin 



2 I Qd 



UUc 

~2J 



< 0(57) 



at the limit of i 



0. 
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II III 




II 



D 



0,5 1 1.5 2 

kd + e 



FIG. 1: Stability phase diagrams for (a) ncU/\J\ = 1, and (b) 
ncU/\J\ = 2 in ID system. The regions S, L and D represent those 
in which the eigenvalue fi^'j^^ is positive, negative and complex, 
respectively. We classify the axis of abscissas into three regions, 
denoted by I, II and III. The eigenvalue is always complex in the 
region II, while it is always real in the region I. In the region III, 
it can be real or complex. 



A. Some properties of eigenvalues and eigenvectors 



We investigate some common properties of the eigen- 
values hujl^^ and eigenvectors a;^^^ . 

We find the foUowing two algebraic properties of the 
non-Hermitian matrix Tkq which may be expressed as 

Tfcq = fn^l^^ao + (Afcq + U nc)CT3 + iU ncCT2 , (58) 

where ctq is the unit matrix and ai {i ~ 1, 2, 3) represents 
i-th Pauli matrix. 

As the first property, the matrix Tkq is pseudo- 
Hermitian: 



Then the inner-product of the doublet eigenvectors is 
naturally introduced as 



(59) 



{x, x') = x'a^x 



(60) 



since we will later see the orthogonality of eigenvectors 
with respect to this inner-product, coming from the re- 
lation 



{x,Tx') = {Tx,x') 



(61) 



for any pseudo-Hermitian matrix T in the sense of 
Eq. . Note that the metric with respect to this inner- 
product is indefinite. 

The second property is derived as 

0-l7feq(Tl = fl^kqO-Q - (Afcq + U nc)(T3 - lU nc<72 

(Afc_q + Unc)a3 - iUnca2 
(62) 

It is also clear from Eqs. ([5^ - ([55)) that the eigenvalues 

(63) 



'■'k-q 
-Tk-q 



^i^^ ^^'^ ^k^q related to each other 



Thus the eigenequations Tk-qX^^}- = f^]k'}qx')^}^ imply 
Tkq<^iX^k}q = ^kq^ix'^it-q^ ^'^'^ therefore 



„(±) 



(64) 



We may introduce the traceless part of Tkq by 

Tkq = (Afcq + U nc)(T3 -I- iUnc(T2 , (65) 
for which we have the eigenequations 

(66) 



^kq-^kq 



where huj^l is defined in Eq. ([5^ . Note that a;^^'' and 
x'^,^- are degenerate states because hcu'^l = f'^k-q- '^^^ 
matrix Tkq satisfies the algebraic relation, 

CTlTkqCFl = -Tkq ■ (67) 

Similarly in the previous paragraph, we easily find that 
Tkqa.xi^^ = T/»45'^i^i? (68) 

leading to 



^kq 



aix 



(±) 

kq 



(69) 



Thus we can express the four eigenvectors ajj,^'' and 

x'"i^- only by two parameters, say Ukq and Vkq, except 
for a trivial phase factor: 



„(+) _ ^(+) _ f "fcq 



™(-) _ ™{-) _ 
•^kq ^ -^k-q ^ 



Vkq 

Vkq 
Ukq 



(70) 
(71) 



B. Real eigenvalues 



In this subsection we consider the case where all eigen- 
values are real, e.g., the region I in Fig. 1. From Eq. (|6ip . 
one obtains the following relation 

(^2^-^15)(-i3>-^^) = 0, (72) 

where s and s' are + or — . Therefore we obtain the 
orthogonality condition 

(4^>4V) = (4t-'4^) = o. (73) 

Noting that Ukq and Vkq are real numbers, we may set 
the normalization condition as 



2 2 

Ukq - V. 



kq 



= 1 



(74) 



Then Eqs. ([70]) and (|7T|) imply that the eigenvectors are 
normalized with the indefinite metric as 



y-^kq ' -^kq ) — ^ ^ 
i^kq^T^kq') = ■ 



(75) 
(76) 
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The elements of the eigenvectors are determined ex- 
pUcitly as 



Ukq 



Vkq 



\ 



(2) 
kq 



\ 




(77) 



(78) 



Note that Uj,q and w^g are singular in the limit of £ — > 0. 
The matrix Tkq is diagonalized by the matrix Okq as 



OkqTkqOkq - 



huj 



(+) 

kq 







huj 



kq 



where 



kq 



™(+) ™(-) 

■^kq -^kq 



(79) 



(80) 



The orthonormal conditions ([75)1 , (1751) ^^^d (|76p are equiv- 
alent to Ol-asOkq = as, yielding 



Ok^='J30y3 

The operators, defined by 

bkq = 



(81) 




^kq ^kq 
^Vkq Ukq 



Ckq 

^k-q 



(82) 



satisfy the bosonic commutation relations [bkq, &^-/ 



(5qq' and [bkq, &fcq'] = 0. By using the 6-operators, one 
obtains the equation of time evolution in a diagonalized 
form from Eq. (|48| . 



iri—bkq[t) = 



(+) 







where Eq. (p5)) has been made use of. Thus the time 
evolution of the operator bkq is simply given as 



bkq{t) = bkqe 



(84) 



In terms of the 6-operators, the unperturbed Hamiltonian 
(Uni) is diagonalized as 



Ho^Y. + const. 



(85) 



The real eigenvalue fu^l^"^ can be interpreted as an energy 
of the quasi-particle, and can be negative in our present 



model. The negative energy causes an energetic insta- 
bility which is called the Landau instability. For one- 
dimensional system, the condition of the quasi-particle 
energy ftw^^"* being negative is written as 



cos^ ( ^ ) > cos(fcd + 6) j cos(fcd + 9) + 



' qd 

sin{kd + 9) sin(^d) < 



"27 



(86) 



at the limit of e ^ 0. 

The operator fc^g represents the zero-mode, as its en- 
ergy eigenvalue becomes 

^fco = (ee)' (2C/nc +ee)' ^ as e , (87) 

from Eqs. (|iS)) . ([5^ - ([55)) . This zero- mode is presumed 
to be the Nambu-Goldstone (NG) mode appearing in the 
spontaneous breakdown of global phase symmetry. In or- 
der to prove that it is actually the NG mode, one needs to 
check the Ward-Takahashi relations [s^ , but this has not 
been confirmed yet when complex modes appear. The 
singular elements of ccj^'' are calculated from Eqs. ((77|) 
and ([75]) as 



UkO 



VkO 



-(4e6-)- 



a + 



Mil 

2a J 

_ Mi 

2a 



+ 0(ei), (88) 



(89) 



where a = (t/nc/2)^^'*. Regularizing the zero-mode this 
way, we can include it as a real mode, as was done in our 
previous works [13, HH . 



C. Complex eigenvalues 

Next, let us consider the case where all eigenvalues 
are complex, e.g., the region II in Fig. 1. Then the 
eigenvalues ff^^"^ are complex conjugate to each other 

as hujj^-' ~ f^kq ■ ^o'^ elements of eigenvector Ukq 
and Vkq become complex numbers. From the relation 

{ft^kJ* ~ ^'^fc*^)(^Cfcq', a^fc^) = where the superscript s 
stands for ±, we obtain 



(+) ™(+) 

kq ' ■^kq 



™(+) 

V-^fcq ' -^kq 



= c. 



(90) 
(91) 



where Eqs. (|7D)) and ([7T|) have been used and C turns 
out to be a pure imaginary constant. In order to fix 
the constant C in Eq. ([9T|) . we take uf,^ — w|q ~ 1 for 
convenience. Then the expressions for the elements of 
the eigenvectors in Eqs. ([77)1 and ([75)1 are true for pure 



(2) 

imaginary fiw)^^, and the simple relation, 



(92) 
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is found. This choice corresponds to fixing C = i, and the 
orthonormal conditions ([50)1 and (pij) are summarized as 

0U<J30k^ = <J2 , (93) 

which derives 

= '^20^3 , (94) 

Similarly as in Eq (|82p . we introduce the new opera- 
tors, 

f Ukq iul- \ / Ckq \ 

(95) 

The operators Akq and Bkq satisfy the following relations 

Akq = ^Al_-, (96) 
Bkq = -iSLq-. (97) 

The commutation relations among the Akq and Bkq op- 
erators become 

= [4,-,^!,-'] = -*v.-'' (99) 

[4.-'4,-']-[^'^.-'^y = 0. (100) 

The equation of time evolution (HS)) is reduced to the 
diagonalized form, 

so the operators Akq and Bj^_- develop in time with the 
complex frequencies as follows: 

Akqit) = e-'^''i''Akq, (102) 



The unperturbed Hamiltonian (j45p is reduced to 



This way the Hamiltonian for complex modes is put into 
a diagonal form, but does not have a representation in 
a Fock space, and the complex eigenvalue can not be 
interpreted as a quasi-particle energy. 

D. Hamiltonian and canonical commutation 
relations 



Generally both real and complex eigenvalues can co- 
exist, e.g., the region III in Fig. 1. In such a case, the 
arguments in the preceding two subsections give the fol- 
lowing Hamiltonian in a diagonal form. 



Bo - 2-^ y—^B^-^Akq^ H A^q^Bkq^j 

+ E W + const., (105) 



where the Bloch wave number q is distinguished by its 
subscript depending on the property of the eigenvalue, 
i.e., Qr for real eigenvalue and Qc for complex one. 

The operator Ckq is now written in terms of bkq, Akq 
and bI_- as 

CkqAt) = Ukq,Akq,e '=<!c -«Ufcq^Bfc_q^e '=lc . 

The expansion of the field operator Eq. (|47p becomes 



I 

H^) = (ukqJkqAx)bkq,e~^'^''^r'' + ^fcq^ /fcqr (^)^I-o.g'"'°~'^'' *) 

+ E ('^'^.-e/fcq-c(a^)^5-ce"*"^'-^* - »<,-Jfefe(a;)^Lq-ee"'"^'-^**) ■ (106) 



One can easily check that the field operators satisfy the CCRs in the first band, 

[${x,t)J^{x',t)] = 5{x-x'), 
[Ux,t),^ix',t)] = [^\x,t),^\x',t)]^0.{107) 
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E. Eigenstate of complex mode 



In order to find eigenstates of the Hamiltonian (|105[) 
we introduce the operator that transforms Ckq to bkq or 
Akq and Bkq, defined by 

l/ = exp[^z0fcq-Gfeq-] , (108) 



where Gkq is 



Gkq — « \ CkqCk-q ^^'kq'^'k-q 



(109) 



Note that Okq is in general a complex number, and it is 
written by a pair of real numbers, and 6\,^, as 



Jkq 



hq 



kq • 



(110) 



so the operator V is not unitary unless the parameter 
^fcq vanishes identically. As is well known, the operator 
V gives rise to the Bogoliubov transformation when all 
the parameters 9kq are real [3^ . Even for complex 9kq, 
the transformation generated by the operator V may be 
written down as 



V 



^kq 
^k-q 



where 



O-H0kq) 



= Q- 



cosh{9kq) 
- sinh(6'fcq) 



ykq) 



^kq 
^k-q 



- smh{9kq) 
coah{9kq) 



(111) 



(112) 



Now, we look for the parameter values of 9kq for which 
the above 0^^{9kq) reproduces the matrix O^^Ckq in 
Eq. ([82| for real eigenvalues or in Eq. ((95|) for complex 
ones. Simple calculations show 





^kq 


= 0, 


(113) 


cosh 


9R 
^kq 


= Ukq, 


(114) 


sinh 


9R 
^kq 


= Vkq 


(115) 



for real eigenvalues, and 

^ (cosh sinh 



^kq 



kq ) = ""fcq 



(116) 
(117) 



for complex ones. 

Let us focus on the complex eigenvalue sector in which 
we have explicitly 



ifcgc 



B 



and 



k — qc 



^kq. 



Bk-q, 



VCkq^V-' 

H-qV- 



^kq, 



i^k—q, 



Ck-q^V^ 



'iBl- 

fcqc 



(118) 
(119) 



(120) 
(121) 



We first define new vacuum states by 

|0)a = T>|0)c, (122) 

|0)b = V-^^\Q)c, (123) 

and 

a(0| = MV\ (124) 

s(0| = MV-\ (125) 

where |0)c is the vacuum of Ckq. They are annihilated by 

^fcQc I ^Iq„ ' -Sfcqc and b],-^ as follows: 

^q-c|0>A = iU-jO)A = 0, (126) 

BkqMB = BlqJ\0)B^O, (127) 



and 



A(0|ifcq-, =A(0|i^ =0, 

B(0|Bfcq-, =b(0|B^-^ =0. 



(128) 
(129) 



From the commutation relations (|55|) - (|100p . it turns out 
that these states are eigenstates of the Hamiltonian Hq, 

,(+)* 

/ X — ^ ' 

M 



H,\0)b = - E 



(+)' 

kqc 



\0)i 



Note the relations for ^(0| and _b(0| , 



^(0|i/o 



B(0|iJo 



E 

V qo 

E 



(+)■ 

kqc 



(+)* 

kqc 



(131) 

(132) 
(133) 



The commutation relations suggest us to intro- 
duce the following excited states through cyclic opera- 
tions of i?t on |0)^ and of A' on \0)b- 

\Nc)a = |?T-q,i • • ■nq^.)A 



Qaj V '^"^ 

\Nc)b = hq^i ■■■nqjB 

^fcqc 



njA(4q-.r-|0)i3. (135) 



These excited states are also eigenstates of Hq as 



(+)* 



ho\nc)a = E "«-c^4£ - 



( + ) ^^fcq, 



I^c)a, 



^0|^c)i3 = E ('^9c^^, 



(136) 



(+)* _ f^kql 

kqc 2 



(137) 



10 



Let us evaluate yi(0|0)A, which is rewritten as 

a{0\0)a = c{0\W\0)c 
W = exp[^|^ 

9o 



(138) 



Consider a state |^), 
\C) ^ exp[^-e 

= exp 1^ - In cosh C + X! ^Uc ^fc-gc ^^'^l^ 1 0) c , 

(139) 

then we have the relation A(0|0)y!i = c(0|^ = — i7r/2). 
Obviously c(0|f) diverges at the limit ^ — i7r/2, so 
does a(0|0)a. Similarly b(0|0)_b is also divergent. 
Using Eqs. (|122p - (|125p . one easily derives 

A(0|0)i3 - 1, (140) 
s(0|0)a - 1. (141) 

The commutation relations in Eqs. ([55)) - (|100p lead to 

A{N'c\Nc)B^X{5nL .n,^., (142) 



b{Nc\Nc)a^X{5^. 



n- rin ■ 



(143) 



Let us rewrite the completeness relation in the complex 
mode sector, using the complete set of \Nc)c, 



Ic = ^ \Nc)c c{Nc 

Nc 



(144) 



where Ic is the identity operator in the complex sector 
and 



\N. 



C)c 



•t \na 

— 1 *^ 
kqcj ' 



|0). 



(145) 



Having V and operate on Eq. (|144p from the left and 
right, respectively, we obtain the completeness relation 
using the states \nq^^ , nq^^jA and b (f^gci i ' ' ' > ^gci I 
as 

Ic - ^V^|A^c)c c{Nc\V-^ 

Nc 

= Y,\Nc)a b{Nc\. (146) 

Nc 

Similarly another relation using the states \Nc)b and 
a{Nc\ follows: 



Ic = ^|A^c)s a{Nc\ 

Nc 



(147) 



One may say that a natural conjugate of |A^c)a is s(-^'^c| 
and vice versa [341 . 



V. PHYSICAL STATES 

In the previous section, we have "diagonalized" the 
unperturbed Hamiltonian including complex eigenvalues 
and have found its eigenstates. The state space is not 
a simple Fock one. We need to impose appropriate con- 
ditions to construct a restricted physical state space. In 
QFT, unstable behaviors of system are described in a sta- 
ble picture such as the Beliaev process. We should now 
establish a stable particle picture specified by the un- 
perturbed Hamiltonian, and the decay processes are de- 
scribed as the higher order of perturbation. We presume 
that unstable behaviors of the BECs in optical lattices 
occur due to external perturbation. 

As in Ref. [s^l , we require the following physical state 
conditions (PSCs). 

i) mi'ix)\n) = v{x), 

ii) (Ti\^\x)^ {x)\Q.) is time-independent , 

iii) (illGjil) is real, when G is any Hermitian operator, 

iv) ^ 1 , 

where (0| is the natural conjugate of jil). If and 
satisfy the above four conditions, we call them physical 
states. The conditions i) and ii) mean that the order 
parameter and density distribution are stationary with- 
out perturbation. The condition iii) guarantees that the 
expectation value of any Hermitian operator can inter- 
preted as physical quantity. The condition iv) is nec- 
essary for the probability interpretation. The vacuum 
states which satisfy the PSCs are obtained as direct sum 
of |0)a and |0)b, 



|0)e ^ -^(|0)a®|0)b) , 
e(0| ^ -L(s(0|©^(0|). 



(148) 
(149) 



The proof that these direct sum states satisfy PSCs is 
given in Ref. [11] . Here we add that the direct sum of the 

excited states lA'^c)^ and \Nc)b are also physical states, 

\Nc)® = -^(|A^c)A©|iVc)s) , (150) 
e(A^c| = -^(B(A^c|©A(iVc|) . (151) 



VI. LINEAR RESPONSE 

So far, we have developed the description of QFT with 
complex eigenvalues. But complex eigenvalues are not 
directly connected with the instability of a condensate. 
In this section, we discuss the dynamics of the system 
with complex eigenvalues, studying the response of a 
condensate against external perturbation. To derive the- 
oretical expressions is straightforward in the linear re- 
sponse theory (LRT) [4§| with our formulation of QFT. 
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We also show numerical results of LRT and compare 
them with those from the TDGP equation, concretely 
those from the discrete nonlinear Schrodinger equation 
(DNSE) which is obtained by applying the tight-binding 
approximation to the TDGP equation [ssllSq. 



A. Formula 

The field operator 5* is expanded in terms of the Wan- 
nier functions as 5" = "^fe + X^i o.kiWki- The particle num- 
ber operator N — jSx'^'^'^ is written as 



i 

Pi{t) - pT{t) + pt\t)^ 



(152) 



(153) 



(154) 
(155) 



where Piit) is the particle density operator at the i-th 
site, and aki{t) is 

1 Ckqr (t) + V 1 Cfcq^ (t) , 

T2 ~ T2 



(156) 



where 



(157) 



„(+). 



Ckq^t) = Ukq.Akq^e -^feg^S^^-^e '='Jc 



(158) 



We consider the external perturbation as 



+ nSe^''-^^al + alakj] , (159) 



where 



SVij{t) = d XwliSVpcrWkj 



(160) 



The function dVpcr{x,t) represents the time-dependent 
modification of trap. We use the on-site approximation 
for SVij {t) . The external perturbative Hamiltonian (|159p 
becomes 



^pcr(i) =^'5yiW/5i(i), 



(161) 



where we write SVi{t) = 5Vii{t). From the linear re- 
sponse theory (LRT) [i^], the change in the i-th site 
particle density {Spi{t)) is given as 

(Spiit)) = {Sp'i'\t)) + (Spt\t)), (162) 
{Spf\t)) = Y.\f dt'{[pf\i),Pi{t')])5Vi{t'), 

j J -CXI 

(163) 

{Spt^t}) = El- /' dt'{[pf-\t),pj{t')])SV,{t'). 

(164) 

Here the expectation (•) is taken to be (A^ | • |iV) where | A^) 
is a direct product of a Fock state for real mode \Nji) and 
the physical state for complex mode |iVc)© in Eq. (I150|) 
with Eqs. (fT34l) and (fTSS]) : 



\N) = \NK)\Nc)e, 
\Nn) = \nq^,---)b 



n 



(165) 
(166) 



The correlation function of Eq. (|163p becomes 



{[pritipjit')]) 

-J^^^iukq+Vkq) 



sin {q-{x,-Xj)-J+\t-t')Y (167) 



Note that q — mode is cancelled in Eq. (|167p and does 
not show divergence in the limit of e ^ 0. 

The correlation function of Eq. p64p becomes 



{[pt\t),Pjit')]) 



2i 



qi,q2 



XRe Ukq^Vkq^{ukq^Vkq^ + Ukq^VkqJ 

X sin|(qi - ^2) ' (xi - Xj) 

This term gives rise to the singularity at q = in the 
limit e — !■ 0. The infrared divergence caused by the zero- 
mode singularity can be removed in the careful treat- 
ments of renormalization [27j or of the quantum coordi- 
nates [2^. In the numerical calculations below we drop 
the divergent term in Eq. (|168p for simplicity, since the 
zero-energy contributions are numerically small after the 
treatments in Refs. [13, Hi. 
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B. Numerical Result 

In this subsection, we show some numerical resuhs of 
LRT with complex eigenvalues and compare those ob- 
tained from DNSE. Here we assume a system of one di- 
mension in space for simplicity. 

DNSE is given as follows 



QOCtE 



at 



(169) 



Here the quantity l^'^P represents the density of conden- 
sate particle at the i-th site. Without external potential 
{5Vi{t) = 0), DNSE has a stationary solution 



(170) 



We adopt this solution for the initial state of the wave- 
function, and define the density response as 



<Sn.(t) = |vI/,(t)|2-|vl/M(<)|2. 
We focus on the external potential of the form 



(171) 



SVi{t) = S-exp 



eit-to), (172) 



which pushes up the center of the lattice, switched on at 
t — tQ. The detailed form of the density response is given 
in Appendix |B] We have calculated the condensate par- 
ticle density n{x) numerically, and have confirmed that 
n{x) is symmetric under the conversion of a; — > — x, i.e., 
all the GP solutions for the one-dimensional system we 
found satisfy n{x) = n{—x). When n{x) has the reflec- 
tion symmetry, the Wannier functions and J become real. 
So we can set 8 = which is the phase of J. 

We set the parameter S/J = 0.0001 , U/J = 0.01 , 
a/d = 0.2 and the total number of lattice sites Is — 51 
with the condensate particle number per site ric = 2. 

The density response in Eq. (|B2p - (|B4[) depends on the 
choice of the state |A^) in Eq. (|165p . As we are inter- 
ested mainly in complex modes in this paper, we take 
the vacuum for real modes in our numerical calculations: 

|iVR,) = |0)b. (173) 

First, we take the vacuum state of complex modes, 

l^c)e = |0)®. (174) 

The critical value of k, giving a boundary between real 
and complex modes, is determined from kd = tt/2. In 
Fig. [51 the time evolution of the changes of the density 
at the 10th site are plotted for (a) kd = IOtt/SI and (b) 
kd ~ 207r/51, in both of the cases all the eigenvalues are 
real. We can see that the results of LRT fit that of DNSE 
with high precision, as is expected. 

Move to the case in which complex modes appear. In 
Fig. [21 the time evolution fo the change of the density 



(a) kd = ^ 




SO 
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T 
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Q0CO1 






-QOOOI 




-QOQOe 
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c 
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-aoue 




-Qoooe 
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-KD 
T 



133 



2D0 



FIG. 2: (Color online) Time evolution for the change of the density 
at the 10th site for (a) kd = 107r/51 and (b) kd = 207r/51 (region 
for real eigenvalues) with T = The solid line represents the 
result of DNSE <5nio(T) while the dashed line represents that of 
LRT (5pio(T)>. 



at the 10th site is plotted for kd = 267r/51, this time 
the complex eigenvalues exist. One can find that the 
density response for small perturbation show the charac- 
teristic behavior. The change of density for kd = 267r/51 
are larger than that of kd = IOtt/SI or kd — 207r/51 
and grows exponentially. This behavior is caused by the 
complex eigenvalues. The result of LRT is in good agree- 
ment with that of DNSE again. Recall that the operators 
^kqc^ Bkq^ and thc state jO)® are essential in our LRT 
formulation. The above agreement is not trivial at all 
when there are complex eigenmodes. 

Next, we focus on one of the complex modes, charac- 
terized by the Bloch vector = 27r/51, and employ the 
following singly excited state for |A^c)e 



\Nc)® = \0---ng^---0), 



(175) 



In Fig. m the quantity {5pio{t)) is plotted as a function 
of time. As riq^ increases, the amplitude of the response 
becomes larger and the exponentially diverging behavior 
becomes prominent. Thus it seems that the excitations 
of the complex modes hasten the collapse of condensates. 
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start with the consistent stable particle in QFT. 

As an application of our formulation to the problem of 
the dynamical instability, we consider the response of a 
condensate against external perturbation in the regime of 
the linear response theory (LRT). The numerical results 
of LRT with complex modes are to be compared with 
those from the discrete Schrodinger equation (DNSE). It 
is remarkable that both of the results coincide with each 
other although the two methods are quite different. The 
state 10)0 and the operators A^g^ and Bfcq^ are crucial for 
our formulation of LRT. It is an interesting observation 
that the excited state in the complex mode sector hastens 
the collapse of the condensate in comparison with the 
vacuum state. 



FIG. 3: (Color online) Time evolution of the change of the den- 
sity at the 10th site for kd = 267r/51 (region for complex eigenval- 
ues) with T = Jt/h. The solid line represents the result of DNSE 
<5nio(T) while the dashed line represents that of LRT (5/3io(T)). 




-Q6 . : I 

30D 3SD 4C0 4SD 5D0 

T 

FIG. 4: (Color online) Time evolution of the density at the 10th 
site ((5/5io(T)) with T = Jt/H, calculated with the state defined in 
Eq. ITTSl l for the interval between T = 300 and T = 500. The solid, 
dashed and dashed-dotted lines correspond to riq^ = 1, 10 and 20, 
respectively. 



VII. SUMMARY 

We have investigated the dynamical instability of Bose- 
Einstein condensation flowing in an optical lattice. The 
formulation in this paper is based on quantum field the- 
ory (QFT). The field operator must include complex 
modes for the dynamically unstable system, otherwise 
the canonical commutation relations would be violated. 

We have attempted to diagonalize the unperturbed 
Hamiltonian under the tight-binding approximation, but 
it is not diagonalizable in the conventional bosonic repre- 
sentation. Nevertheless one can find its eigenstates, i.e., 
the vacuum and excited states in the complex mode sec- 
tors. Then, the physical state conditions (PSCs) were 
introduced to restrict the state space, so that we can 
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APPENDIX A: EIGENFUNCTIONS OF THE 
BOGOLIUBOV-DE GENNES EQUATION 

In this Appendix, we rephrase the contents in Sec. IIVI 
from the viewpoint of the Bogoliubov-de Gennes (BdG) 
equation. The QFT formalism on the trapped BECs with 
vortices in Ref. [SJ] is based on the BdG equation. 

The relevant BdG equation is, in the doublet notation, 

TXn{x) = hoJnXnix) (Al) 

where 

.4-) = (A3) 

L = K + Vopt- f^ + 2g\v{x)W (A4) 
M = gv'^ix). (A5) 

The operator T has the pseudo-Hermitian property of 

asTtaa = T (A6) 

correspondingly to Eq. ([SS]) . This leads us to define the 
following inner product for an arbitrary pair of doublets, 

(r, s) = J tfxr''{x)a3s{x) (A7) 
= Jd^x[rl{x)si{x)-r;{x)s2{x)], (A8) 
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where 



where 



r{x) 



ri{x) 
r2ix) 



six) 



si{x) 



(A9) 



We may also define a (squared) "norm" of r as = 
(r, r), which is not positive-definite. One easily obtains 

{hWn'-hw;^,)iXn',Xn) ^0. (AlO) 

As a counterpart of Eq. ((62l) we find the relation, 

aiT*ai = -T. (All) 

It turns out that for any eigenvector Xn{x) whose eigen- 
value is denoted by hu>n the doublet yn{x) ~ aix*-^{x) 
becomes an eigenvector with the eigenvalue —fujJn- 

When the eigenvalues are real, we have the following 
orthonormal relations, consistent with Eq. (jAlOP : 



(Un'^Xn) = 0. 



(A12) 
(A13) 
(A14) 



Complex modes appear in a pair for the BdG equa- 
tion (|A1[) . i.e., any eigenstate belonging to a complex 
eigenvalue hcom is accompanied by another eigenstate 77^ 
whose eigenvalue Huji is a complex conjugate of fkOm, 
fkoi = hu^. This fact is shown in constructing eigen- 
states explicitly in Sec. IIVI The "norm" of the eigen- 
states of complex eigenvalues is zero. 



lie™ II 



= 0. 



(A15) 



since {hujm — zero in Eq. (jAlOp . The "zero 

norm" is a necessary condition for the emergence of com- 
plex eigenvalues. The pair of the eigenvectors Cm and rfi 
are not orthogonal to each other in general. 



ive,Cm) ^ 



(A16) 

huje) on the left- 



as there is a vanishing factor of (fiw^ 
hand side of Eq. (|A10|1 . 

Let us expand the field operator in the doublet nota- 
tion in terms of the eigenfunctions of the BdG equation 

[Cfcq. {x)Akq^ (t) + r]kq^ (x) BI-^ (t) 



E 



(A17) 



Hx,t) 



( kx,t) 



(A18) 



\4'Hx,t) 

as in the case of a BEG with a vortex [34|. Comparing 
this expansion with Eq. (jl06p and its Hcrmitian conju- 
gate, we have 



Xhqr {x) 



Ukqjkq,{x) 
VkqJk-nAx) J ' 



11, - (■r) — i '"kqrfkq^ix) 
y-"'^''' - I U^qJk^,Sx) 

Ukq^fkqSx) 
-^K-qJl-qS^) 



CkqSx) 



I UkqM,S^) 



(A19) 

(A20) 
(A21) 
(A22) 



It is straightforward to check the following orthonormal 
relations, 



and 



(ykq'^iVkqr) = -^qrq'„' 
{Vkq'^^^kq,) = 0, 



iCkq'^XkqJ = 0, 
{Vkq'^,VkqJ = 0, 
iVkq'^XkqJ = iSq^q'^ 



(A23) 

(A24) 
(A25) 



(A26) 
(A27) 
(A28) 



APPENDIX B: EXPRESSION OF DENSITY 
RESPONSE 

In this Appendix, we give the detailed expression of 
the density response {6pi{t)). The state by which the 
expectation (•) is taken is found in Eq. (|165p . We restrict 
ourselves to the case of one dimension in space. 

The external perturbation is given as 



SV^{t) = S-exp 



Xi - Isd/2 



e{t~h). 



(Bl) 



Then the expression of the density response becomes 

{5Mt))^{5pf\t)) + {5pt\t)), (B2) 

where 
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2S V- 



Xj - Isd/2\ ndukq + V,,q) 



u> 



,(+) 

kq 



COS 



- Xj) - ojl^\t - to)} - cos{q{xi - Xj)} 



(B3) 



T^Y. exp(^^— {(nq, +l)(n_g^ + l)-ng,n_,-J 



hi 

xRe 



j gi ,927^0 

Ukqi Vkq2 (Ukqi Vfeqa + Ukq2 ^kqi ) 



( + ) + ) U 



{ COS ((gi - g2)(ari - a;,) - {J+l + 4-k)(* " *o)) 



cos {{qi - q2){xi - 



(B4) 
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